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1.0 Introduction 

In fitting data, some assumptions about the process to be modelled are in¬ 
evitably made, either implicitly or explicitly. There is much empirical evi¬ 
dence that many naturally occurring processes are modelled well by fractional 
Brownian Motion (fBm), and this model can form the basis for a model-based 
statistical estimator, the fBm interpolator. 

While interpolation is often considered for the case of doubles {x,z) or 
triples (a:,?/, 2 ), the derivation that follows is easily generalized to data of 
higher dimension. For this reason, it is convenient to represent position as a 
vector (e.g., a = {x,y)). 

Given n data points of the form (a,, z,)? the interpolation problem can be 
considered as an estimation problem: find the best estimate of z at the location 
a. Obviously, finding the best estimator using classical methods requires a 
detailed knowledge of the statistics of the process. 


2.0 Statistical Model of the Data 

2.1 Error-Free Data 

Let us assume (a,-,z,), z = 1,..., n to be n error-free samples of a fractional 
Brownian motion. While we assume no error in the data (i.e., the data values 
are exact), this assumption is probably not realistic, and it is relaxed in section 
2 . 2 . 
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If (aijZi) is fBm, then it is well known that the increments Zi — Zj are 
Gaussian with zero mean and variance proportional to the spacing of the 
points: 

va.T{zi - Zj) - a'^Wai - ( 1 ) 

where || ■ || represents the usual Euclidean norm. We assume the statistics to 
be locally homogeneous, that is H and are assumed locally constant. The 
exponent H describes the roughness of the surface, and it must be estimated. 
It takes values in the range 0 < H < 1, where H = 0.5 corresponds to true 
Brownian motion, and empirical evidence indicates that 0.7 < H < 0.99 corre¬ 
sponds to many natural phenomena. The scale factor can be interpreted as 
the mean-squared “variability” of the fBm. More detail on the mathematical 
foundations of fBm can be found in [1]. 

To simplify the notation that follows, we define the quantities 

Pi - ||ai - and pij = Ija^ - ■ (2) 

Thus we can write var(zi — z) = c'^pi and var(z,- — Zj) = Furthermore, 

since z is fBm, 

cov(zi,Zj) = i<r^l2)(pi + Pj - pi,j). (3) 

Let z be a vector of data values, z = [zi] for i = 1,..., n. Since z and the 
interpolation point (a, z) are sampled fBm, the data vector z — zl (where 1 is 
a vector of all ones) is multivariate Gaussian with zero mean and covariance S. 
Viewed another way, z is multivariate Gaussian with mean zl and covariance 

S. 

Assuming the unknown z value is fixed, the covariance matrix is given by 

S = Ezz'-z^ll' (4) 

= cr^R 

where the matrix R has elements = {pi -f pj — pi^j)/2. Note that the 
covariance is factored into R a function of data and interpolant positions, and 
cr^. In the error-free case, we can generate an estimator that is independent of 


2.2 Data with Additive Errors 

Assume the data z to be a fBm vector z with additive Gaussian error e, that 
is z = z -f e. The covariance of z can be easily found. The statistics of the 
fBm vector z are given in Equation (4) for the error-free case. Let the error 
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have zero mean and covariance Eee' = rj^C and independent of the fBm, z. It 
follows that the covariance matrix of z is given by 

S = cr^R -f if'C (5) 

where R is defined in Equation (4), and C describes the mean-squared error 
in the data. 

In addition to cr^, this introduces a new quantity, the error covariance, r/^C. 
Unfortunately, unlike the error-free case, we cannot generate an estimator 
that is independent of 77 ^, cr^, and C. Fortunately, it is frequently possible to 
determine the covariance structure. 

Several cases come immediately to mind. In the simplest case, we assume 
iid measurement errors with zero mean and variance 7 ;^, then we can write 
C = I, and E = cr^R -f t;^!. The resulting estimator requires one additional 
parameter not needed in the error-free case: 

While this may be true of some data sets, frequently, as in the case of 
multibeam sonar bathymetry, there are systematic and unavoidable errors in 
the data collection that produce different levels of error at different points. In 
this case, C might be assumed to be diagonal {independent errors at adjacent 
points), but not constant on the diagonal. Merging data sets with different 
mean-squared errors will produce the same result. Finally, in the most general 
case C need not be diagonal. 


3.0 The Mciximum Likelihood Estimator 


Viewed from the position a, fBm can be expressed as a Gaussian vector with 
constant mean zl and covariance S. Thus the estimation of z at a can be re¬ 
stated as the estimation of the mean of a Gaussian vector. This is a well known 
problem in statistics and the Maximum Likelihood Estimator is presented here 
without proof. For further information see Mardia et al.[2]. 

Given a multivariate Gaussian with mean zl and covariance S, the maxi¬ 
mum likelihood estimate of z is easily found to be [ 2 ] 


rs-^z 

I'S-il 


w'z. 


( 6 ) 


The estimator z is an unbiased minimum-variance estimator with expected 
value 2 : and variance 

var(z)/cr^ = 1/1'S“U. (7) 

Note that the weight vector w is a function of data positions 5, and not of 
z-values. 
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3.1 Error-Free Data 


For the special case of error-free data, the covariance of the data vector z was 
shown above to be S = cr^R. Thus, the maximum likelihood estimator of 
Equation (6) becomes 

■ ■ ( 8 ) 


I'R-U 


= w z 


As promised, the estimator does not require knowledge of cr^. The normalized 
variance ^ of the estimator is 


= var(i)/cr^ = 1/l^R ^1 


(9) 


The normalized variance will prove to be a useful measure of performance that 
is independent of cr^. 


3.2 Data with Additive Errors 

For data with additive measurement error, the maximum likelihood estimate 
of 2 in Equation (6) becomes 

. T(R + i/2c)-i2 

r(R-fi/2C)-U 

where The estimator is unbiased with expected value 2 and 

minimum normalized variance 

^ = var(2)/(j^ = l/l'(R + (11) 

Note that the weight vector w is a function of data positions a,- the ratio 
77 ^/cr^, and the error correlation C, but not of 2 - values. 



4.0 Iterative Estimation of z 

The estimator given in (6) can be implemented iteratively. This follows imme¬ 
diately from the observation that the n-point estimator 2 („) and the (n -f l)-st 
data point are jointly Gaussian. While this may not reduce the computational 
burden of this routine given the high efficiency of the Cholesky decomposition 
that can be used in the direct solution of Equations (8) or (10), there are two 
potential benefits of this technique. Since the variance of the estimator is com¬ 
puted at each stage in the algorithm, a stopping rule can be used to terminate 
the iteration after a predetermined level of performance is achieved. At each 
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stage, multiple points can be tested, the resulting performance evaluated, and 
that point which most decreases the variance of the estimator can be selected. 
Using the iterative algorithm combined with these rules produces a low order 
estimator based on a small, carefully selected set of points with a reasonable 
level of performance. 


4.1 Derivation of the Iterative Algorithm 


Let Z(n) — be the estimate of 2 based on the length n data vector 

Z(„). As stated above, it is Gaussian with mean z and variance An 

updated estimate of z can be found using Z(„) and an additional data point 
((3ri,+i,-2n+i)- A new data vector can be made from the n-th order estimate 
and a new data point; this yields the vector [z(„) z„+i]' with mean [z z]' 

and covariance 


^(n+l) — 


({n) 

L ^(n)W(n) 


^(n)W(„) 

Pn+H 


( 12 ) 


where 


r(„) =cou(z„+i,Z(„)) = 0.5[p,f = l,...,n. (13) 

As above, the MLE of z is given by 


-^(n+l) 


I'E-^z 

rs-u 


1 

Pn+l *'(n)^(n) 

/ 

^(n) 

7 

r 

vr 

1 

_ 1 




where 

7 = ^(n) + p„+i -2r;„)W(„). 

The new normalized variance is 


(14) 

(15) 


'^(n+i) = det(R )/7 

where 

det(R) = ^(„)p„+i - (r(„)W(„))A 
Finally, the weights can be updated as follows: 

1 

W(„+i) — — 

7 




(16) 

(17) 


(18) 
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4.2 The Iterative Algorithm 

The iterative algorithm is presented below. 

0. Assume the scalars and ^(n)) S’Hd the vectors and pi for i = 
1,..., n are known. 

1. Get the new data point (5„+i,iTn+i)- 

2. Compute pn+i and pi^^+i for i = 1,... ,n using Equation (2). 

3. Compute r(„) using Equation (13). 

4. Compute ^(n+i) using Equation (16). 

5. Compute Z(n+i) using Equation (14). 

6. Update the weights W(„ 4 .i) using Equation (18). 

7. Go to 1. 

5.0 Example 

As a simple example of the character of the fBm interpolator, consider a simple 
one dimensional interpolation with data points (x,z) = {(0.2,8.0), (1.6,6.0), 
(2.9,8.0), (4.2,11.0), (5.7,13.0), (7.0,13.5), (8.1,14.5), (9.6,15.0)}. 

The results of interpolating this data are shown in Figures 1-5. The 8 data 
points are interpolated on the grid x = 0 : 0.1 : 10 using different values for H 
and The data points are indicated in each plot by the asterisks. 

In Figures 1 and 2, the data is assumed error-free, thus the interpolating 
function passes through the data points. Note however the markedly dilferent 
behavior of the interpolating function as H is varied. For H = 0.3 in Figure 
1, the function is “peaky,” with cusps at the data points, returning rapidly 
to a mean value, and containing significant high frequency components. This 
corresponds to interpolating a surface with a high fractal dimension (and high 
dimension results in what we instinctively call a noisier surface). Most natu¬ 
rally occurring surfaces or functions tend to have much lower fractal dimension: 
H = 0.9 is a more typical value. In Figure 2, H — 0.9, and the result is a 
smoother function with dominant low frequency components. 

Contrasting Figures 2 and 3 shows the effect of assuming additive mea¬ 
surement error. In Figure 2, the data is assumed error-free, in Figure 3 it is 
assumed to have variance 1.0. Note that the interpolating function in Figure 3 
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does not pass directly through the data points, and it has significantly less cur¬ 
vature. Furthermore, as is increased, the interpolating function approaches 
a straight line. 

In Figure 4, i/'^ is varied from point to point. The point (7.0,13.5) is assigned 
a variance 10 times that of the other points. As a result, the data point in not 
as tightly fit as the other points. 

In Figure 5, the normalized variance of the estimator var{z}/i/^ is plotted 
versus the x values for the interpolation shown in Figure 4. As expected, the 
performance decreases as the interpolation point moves farther from the data 
values. Note that the variance never goes to zero due to the measurement 
error, and that it increases at both ends of the data. It is interesting to note 
the variance in the vicinity of the “bad” point (7.0,13.5). While the variance 
does increase in this region, it is still relatively small due to the contibutions 
of “better” neighboring points. 

Several conclusions can be drawn from this example. First, the exponent 
H must be estimated, assumed, or derived for a reasonable fit to the data. 
Techniques for estimating H or equivalently fractal dimension are still an open 
question. The results of a fit using a nonzero measurement error will very not 
only reduce the mean-squared error of the estimator, but because it tends 
to “straighten” the interpolating function, it should reduce the severity of 
artifacts resulting from exact fits. Finally, through the use of different values 
of at each point, the fit can be selectively “relaxed” at suspected or known 
bad points. 

6.0 Summary 

The maximum likelihood interpolator for fBm is presented above for data 
satisfying this well-defined statistical model. While this technique has its own 
set of problems, they are quite different from those of many interpolators 
currently in use. 

Setting the measurement error equal to zero results in an interpolating 
function that passes through all data points, but the character of the fitting 
function can be modified by selection of the exponent H. This leads to an 
underlying problem, the correct selection of H. Presumably, this selection or 
estimation problem can best be understood in the context of a particular fBm 
process, and techniques for estimating H must be developed in that context. 
As an aside, choice of H affects the spectral character of the resulting interpo¬ 
lating function: small H results in a “high-pass” fit, large if in a “low-pass” 
fit. 
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Assuming nonzero additive error relaxes the fit as expected. This also has 
the effect of reducing the curvature in the fit. Assigning different mean-squared 
error values to data points allows data of different precision to be merged, and 
selectively relaxes the fit at poor data points. The value of mean-squared error 
at the data points is obviously data specific, and presumably it is known or 
can be approximated. In absence of this knowledge, a global variance can be 
used. 

Although the iterative estimator derived above may not may not reduce 
the degree of computation required, it has two significant benefits. First, a 
stopping rule can be used to terminate the iteration after a predetermined level 
of performance is achieved by using the variance of the estimator computed at 
each stage in the algorithm. Second, the data point which most reduces the 
variance of the estimator can be selected at each stage of the iteration. 
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Figure 2: fBm interpolation {H = 0.9, — 0.0) 


10 






12 





